Spectral theory of metastability and extinction in birth-death systems 
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We suggest a general spectral method for calculating statistics of multi-step birth-death processes 
and chemical reactions of the type mA — > nA (m and n are positive integers) which possess an 
absorbing state. The method employs the generating-function formalism in conjunction with the 
Sturm-Liouville theory of linear differential operators. It yields accurate results for the extinction 
statistics and for the quasi-stationary probability distribution, including large deviations, of the 
metastable state. The power of the method is demonstrated on the example of binary annihilation 
and triple branching 2 A — >0, A — > 3 A, representative of the rather general class of dissociation- 
recombination reactions. 

PACS numbers: 05.40.-a, 02.50.Ey, 82.20.-w, 87.23.Cc 
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Since the pioneering works of Delbriick 0, Bartholo- 
may [2j and McQuarrie -3| , kinetics of systems of birth- 
death type, containing a large but finite number of agents 
(such as molecules, bacteria, cells, animals or even hu- 
mans), have attracted much attention in different areas 
of science and become a paradigm of theory of stochastic 
processes 0, . Birth-death models are extensively dis- 
cussed in chemistry, astrochemistry, epidemiology, popu- 
lation biology, cell biochemistry, etc. They are also well 
known in non-equilibrium physics, and can be viewed in 
the context of reaction-limited kinetics on a lattice, as 
opposed to more extensively studied diffusion-limited ki- 
netics |6(. While the behavior of the average number 
of particles in such systems may be describable, at not 
too long times, by (mean-field) rate equations, fluctua- 
tions may lead to important quantitative or even quali- 
tative differences. This necessitates using the more gen- 
eral master equation which deals with the probability of 
having a certain number of particles of each type at time 
t. The master equation is rarely solvable analytically, 
and various approximations, often uncontrolled, are in 
use 0, 0| , such as the Fokker-Planck equation which may 
suffice unless one has to deal with large deviations or ex- 
tinction phenomena 0, @, ■ Not much is known beyond 
the Fokker-Planck description, though in particular cases 
the statistics were determined b y using appro ximations 
in the master equation [j H(l ITlL Till Im Hi Jl^ or, alter- 
natively, by introducing a generating function 0, 0, 0] , 
see below, and developing different approximations in the 
equation describing its evolution [j Uld.!!^ ]. 

In this work we advance the generating function tech- 
nique by marrying it with the Sturm-Liouville theory of 
linear differential operators. This yields a general and 
robust spectral formalism, capable of providing accurate, 
and often analytical, results for extreme statistics in a 
variety of (not necessarily single-step) birth-death sys- 
tems and chemical reactions. We demonstrate the power 
of our method by a simple reaction of binary annihila- 
tion and triple branching. An example of such a reaction 
is recombination of two atoms A and dissociation of the 
molecule A^. A + A — > A^, and A2 + A — ► 3 A, assum- 



ing that the A2 molecules are always at hand ^^|- For 
H or N atoms this reaction occurs at high temperatures 
. We calculate the extinction probability as a function 
of time, the mean time to extinction and the complete 
quasi-stationary probability distribution of the long-lived 
metastable state, intrinsic to this problem. 

Rate equation, master equation, generating function 
and spectral theory. Consider the binary annihilation 
and triple branching reactions 2A -A 0, and A -A 3A 
where \x, A > are rate constants. The rate equa- 
tion (or the mean field theory) of this simple system, 
dn/dt = 2 An — {in 2 , describes a nontrivial attracting 
steady state n s = 2f2, where fl = A//i. Fluctuations, 
caused by discreteness of particles, invalidate this mean- 
field result owing to the existence of the absorbing state 
n = 0: a state from which there is a zero probability 
of exiting. At ^> 1, however, a long-lived (and there- 
fore quasi-stationary) fluctuating metastable state is ob- 
served, once the initial population is not too sparse. The 
statistics of the quasi-stationary state and of the extinc- 
tion times are the subjects of our interest here. 

To account for discreteness of particles, we assume that 
the evolution of the probability P n (t) to find n particles 
at time t is described, for n > 1, by the master equation 

d ;Pn(t) = £[(n + 2)(n+l)P n+2 (t)-n(n-l)P n (f)] 



dt 



+ A[(n-2)P n _ 2 (t)-nP n (*)] 



Let us introduce the generating function 



G(x,t) = Y j x n P n {t), 



(1) 



(2) 



where x is an auxiliary variable. G(x, t) encodes all the 
probabilities: 

1 d n G(x,t) 



P n (t) = 



(3) 



Obviously, G(x = l,t) = 1. Equations and J2J) yield 
a partial differential equation (PDE) for G(x,t): 



dG fi.^ 2 ,d 2 G . . o ,,dG 



(4) 
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As the reaction we are dealing with conserves parity, 
G(x, t) can be written as 

G(x, t) = Ci_G even (x, t) + C 2 G odd (x, t) , (5) 

where ci — *Y^o PinJ^ = 0)> an( i c 2 — 1 — ci- Therefore, 
Geven{x = ±l,t) = 1 and G od d{x = ±l,t) = ±1. The 
steady state solution of Eq. 10} is 



G st (x,t) =ci + c 2 



er/i (v^s) 



erfiiVQ,) 



(0) 



where erfi(x) — (2/y / 7r) e* (ft. Let the number of par- 
ticles at i = be even: n = 2fc , where ko is integer. 
In this case the parity conservation yields c\ = 1 and 
c 2 = 0, so G s t(x) = 1, and the only true steady state is 
the empty state: Po = 1, while the rest of P n are zero 

To see how the population of no = 2ko particles at 
t = becomes extinct, we introduce a new function 
g(x,t) — G(x,t) — G st (x) — G(x,t) — 1 which obeys 
Eq. (@J with homogenous boundary conditions g(x = 
±l,i) = 0. Substituting g(x,t) = e~ lt <p(x), we obtain 

(1 - x 2 )ip"(x) + 2Qx(x 2 - 1 VO) + 2Etp(x) = , (7) 

where E = j/fj,. Rewriting this ordinary differential 
equation in a self-adjoint form, 



[<p'[x) exp(-ttx 2 )] ' + Ew(x)<p(x) = , 



(8) 



with the weight function w(x) — 2 exp(— Q,x 2 )(l — x 2 )^ 1 , 
we arrive at a standard eigenvalue problem of the Sturm- 
Liouville theory |2l|. Once we have found the com- 
plete set of orthogonal eigenfunctions tpk(%) (which are 
all even), and the real eigenvalues Ek, k = 1,2,..., we 
can solve the time-dependent problem: 



G(x, t) = 1 + } a k (p k (x)e 



k=l 



(9) 



where 



J^[G(x,t = 0) — l]<fk(x)w(x)dx 
Jo ft( x ) w ( x ) dx 



(10) 



,.2A- ( , 



and G(x, t = Q) = x" 

As all Ek > 0, Eq. JHJ describes decay of initially pop- 
ulated states k = 1, 2, . . ., and the system approaches the 
empty state G(x,t — > oo) = 1. We are interested in the 
case of f2 ^> 1, where the metastable state is expected to 
be long-lived. Elgart and Kamenev showed that the 
eigenvalues E 2 , E 3} . . . scale like 0(Q) ^> 1. In contrast 
to these, the "ground state" eigenvalue E\ is exponen- 
tially small, as will be proved a posteriori. We will be 
interested in sufficiently long times fiVLt = Xt 3> 1, when 



the contribution from the "excited" states to G(x, t) be- 
comes negligible, and we can write 



G(x,t) = l + ai<pi(x)e-» Elt . 



(11) 



Let us proceed to the ground state calculations. 

Ground state calculations. As <pi{x) = <p(x) is an even 
function, it suffices to consider the interval < x < 1. 
We will employ the strong inequality 1]>1 and find the 
(very small) eigenvalue E\ and the corresponding eigen- 
function of Eq. Q by a matched asymptotic expansion, 
see e.g. Ref. |2jj. In most of the region < x < 1 (the 
bulk) we can treat the last term in Eq. Q perturba- 
tively. In the zero order we put E\ = and arrive at the 
steady state equation tp' b '{x) — 2VLxip' b {x) = 0, whose even 
solution is <p^ — 1. Now we put ifb(x) = 1 + 5(fb(x), 
where dipb(x) <C 1, and obtain in the first order 

5ip%(x) - 2Q.x6<p b (x) = -2J5 X (1 - x 2 )- 1 . (12) 

The solution for ifib{x) takes the form: 

f x f s e~ 

<p b (x) = l-2E 1 / e ns ds / ^ dr ■ ( 13 ) 

Jo Jo 1 — r 

Asfi > 1 , we can omit the r 2 term in the denominator of 
the inner integral in Eq. (this omission is illegitimate 
too close to x = 1, but the bulk solution is invalid there 
anyway, see below). We obtain 



(f b (x) ~ 1 - 2Ei \ e ns ds / e-° r dr 
Jo Jo 

= l-E lX 2 2 F 2 (l,l;^,2;ihA , (14) 



where 2 F 2 (ai, a 2 ;bi,b 2 ;x) is the generalized hypergeo- 
metric function |2J] , while Ei is yet unknown. It is easy 
to check that the bulk solution is valid [S(fb (x) <C 1] as 
long as 1 — x ^S> I /ft. 

In the boundary layer 1— x <C 1 we can again disregard, 
at ^> 1, the (exponentially small) last term in Eq. J7|). 
The resulting equation is again <p'/(x) — 2Qxtp' l (x) = 0. 
Its non-trivial solution, obeying the boundary condition 
at x = 1, is 



ipi (x) = C 



ds ^ C 2n 



,-0(l-x 2 ) 



(15) 



where C is a yet unknown constant. 

Now we demand that, in the common region l/Q, <C 
1 — x <C 1, the proper asymptote of the bulk solution 
(|14f) . obtained by moving to infinity the upper limit in 
the inner integral of Eq. (|14|) : 



V^Ei 



ds ~ 1 



2r>3/2 



(16) 
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coincides with the boundary layer solution JTSJ. Equa- 
tions iJTSJ and (jTBjl yield 



and C = 2fte~ 



(17) 



As expected, the lowest eigenvalue E\ is exponentially 
small in ft. The respective eigenfunction is 



ip{x) ~ < 



' <p b (x) = I - E lX 2 2 F 2 (1,1; |,2;fta; 2 ) 
for 1 - x 2 > 1/ft, 

pi (a:) = l- e -"d-^ 2 ) 
for 1 - x 2 < 1 . 



(18) 



Now we use Eq. (|1U|I to calculate the coefheient a± enter- 
ing Eq. iJUJ. While evaluating the integrals, we notice 
that the main contributions come from the bulk region 
1 — x 2 ^> 1/ft, and it suffices to take the eigenfunction 
tpb(x) in the zeroth order, that is ip?\x) ~ 1. Evaluat- 
ing the integral in the nominator of Eq. I|1U|I . we notice 
that, the term x 2k ° under the integral is negligible com- 
pared to 1. As a result, the integrals in the nominator 
and denominator become identical up to a minus sign. 
Therefore, a\ ~ — 1 which completes our solution 1)110. 

Statistics of the quasi- stationary state. We start this 
part with calculating the average number of particles 
n and the standard deviation a at intermediate times 
ft -1 < fit < E-f 1 . Using Eq. {TTJ, we obtain 

n = d x G\ x=1 = 2ft , (19) 

which coincides with the mean field result. Furthermore, 



2 2 -2 

a = n z — n 



= [d 2 xx G + d x G-(d x G) 2 ]\ x=1 = m. 



(20) 

where we have used for tp(x) its boundary layer asymp- 
tote fi(x) from Eq. (|18fl . One can see that, at inter- 
mediate times ft^ 1 <C /it <C E^ , the system stays in 
the (weakly fluctuating) quasi-stationary state. What is 
the complete probability distribution P n (t) of the quasi- 
stationary state at these times? For n — we obtain 



P Q (t) = G(x = 0,t) = l-e~ flElt 



(21) 



which, at \iE\t <C 1, is very small. For (even) nonzero 
values of n, Eqs. © and l(TT)> yield 



, , 2^ 1 (4ft) n / 2 " 1 (n/2- 1)! lFI 
P n {t) = — - — - - e~^ Elt . 



(22) 



For n»l we can use Stirling's formula and obtain 



Pn(t) 



\/2nft 



f(l+ln^)-/iBi« 



(23) 



Notably, all of the probabilities P n {t) (n > 0) decay 
with time, while Po(t) grows. One can check that the 



most probable state coincides with n = 2ft. In the vicin- 
ity of n = n, P n (t) from Eq. (|23|l can be approximated 
by a normal distribution with the mean n and standard 
deviation a, given by Eqs. ifHfl) and (|2*U|> . respectively. 
The tails of the true distribution, however, are strongly 
non-Gaussian. A comparison between our analytic re- 
sult (|22|l . the large-n approximation l|23() . and the nor- 
mal distribution is shown in Fig. ^ One can see that 
Eq. H23J) is very accurate, whereas the gaussian approxi- 
mation strongly overpopulates the low-ro tail and under- 
populates the high-n tail. 
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FIG. 1: (Color online) The natural logarithms of the prob- 
ability distribution 1221 (the blue solid line), of its large-n 
asymptotics (ll'Mt (the green dotted line), and of the normal 
distribution with n and a from Eqs. 119H and 1201 (the red 
dashed line), for Q = 40 and fiEit <C 1. 

Figure [21 compares our analytic result l|22|l with a nu- 
merical solution of the (truncated) master equation 
with (d/dt)P n (t) replaced by zeros and Pq = 0. The two 
curves are almost indistinguishable for ft = 10. In fact, 
good agreement is observed already for ft = 0(1), and it 
rapidly improves further as ft increases. 




FIG. 2: (Color online) The red dashed line: the natural 
logarithm of the probability distribution 1221 at fiEit <C 1. 
The blue solid line: a numerical solution of the master equa- 
tion Q, see text. The parameters are fi = 10, no = 2ko = 40. 



4 



Extinction time statistics. The quantity Po(t) is the 
probability of extinction at time t. The extinction prob- 
ability density is p(t) = dP (t)/dt. Using Eq. J2U, we 
obtain an exponential distribution: 

p{t) ~ fiE ie -^ Elt at Xt > 1 . (24) 

The average time to extinction, f = J °° tp{t) dt ~ 
(/xi^i ) — 1 , is exponentially large, at SI 3> 1, as expected. 
Figure [3] compares the analytical result 112 1H for Pq (i) 
with G(0, i) found by solving Eq. |0J numerically with 
the boundary conditions G(±l,i) = 1 and the initial 
condition G(x,t = 0) = x 2fc °. The inset compares the 
analytical and numerical ground state eigenvalues, and 
good agreement is observed. 



systems of this type the rate equation yields a stable non- 
empty steady state, but an account of fluctuations brings 
about an unlimited population growth, see e.g. Ref. 
In that case the ground-state eigenvalue is expected to 
be negative. Finally, the use of spectral formalism is not 
at all limited to systems possessing an absorbing state 

0. 
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FIG. 3: (Color online) Extinction probability Po(t) from 
Eq. 1211 (the blue dashed line) and from a numerical solution 
of Eq. @ (the red solid line) for Q = 25 and n = 2k = 100. 
The inset shows, at times [it l/fi, the ratio of the numeri- 
cal ground state eigenvalue E™ um = — log[l — G(0, t)]//j, and 
the analytical one, given by Eq. 1171 . The deviation is less 
than 4%, that is within error O(l/0,). 

Final comments. The spectral formalism yields ac- 
curate extinction time statistics and a complete quasi- 
stationary probability distribution of the metastable 
state for a wide class of birth-death processes which pos- 
sess an absorbing state and are describable by a mas- 
ter equation. In this formalism, the problem of com- 
puting these statistics is reduced to a problem (familiar 
to every physicist) of finding a ground-state eigenfunc- 
tion and eigenvalue of a linear differential operator. We 
have demonstrated the formalism by an example of bi- 
nary annihilation and triple branching, but the formalism 
is general and can be used for a variety of kinetics. In 
most interesting cases of long-lived metastable states, a 
large parameter (the average number of particles in the 
metastable state) is always present in the problem. This 
paves the way to a perturbative treatment, like in the 
example we have considered. 

The spectral formalism should be also efficient when 
the absorbing state is at infinity, rather than at zero. For 



[1] M. Delbriick, J. Chem. Phys. 8, 120 (1940). 

[2] A.F. Bartholomay, Bull. Math. Biophys. 20, 175 (1958). 

[3] D.A. McQuarrie, J. Appl. Prob. 4, 413 (1967). 

[4] C.W. Gardiner, Handbook of Stochastic Methods 

(Springer, Berlin, 2004). 
[5] N.G. van Kampen, Stochastic Processes in Physics and 

Chemistry (North-Holland, Amsterdam, 2001). 
[6] J.L. Cardy and U.C. Tauber, J. Stat. Phys. 90 (1-2), 1 

(1998); U.C. Tauber, M. Howard and B.P. Vollmayr-Lee, 

J. Phys. A 38, R79 (2005). 
[7] B. Gaveau, M. Moreau, and J. Toth, Lett. Math. Phys. 

37, 285 (1996). 
[8] V. Elgart and A. Kamenev, Phys. Rev. E 70, 41106 

(2004). 

[9] C.R. Doering, K.V. Sargsyan, and L.M. Sander, Multi- 
scale Model, and Simul. 3, 283 (2005), and references 
therein. 

[10] C.A. Brau, J. Chem. Phys. 47, 1153 (1967). 

[11] W. Nadler and K. Schulten, J. Chem. Phys. 82, 151 

(1985); Z. Phys. B 59, 53 (1985). 
[12] M.I. Dykman, E. Mori, J. Ross, and P.M. Hunt, J. Chem. 

Phys. 100, 5735 (1994). 
[13] I.J. Laurenzi, J. Chem. Phys. 113, 3315 (2000). 
[14] I. Nasell, J. Theor. Biol. 211, 11 (2001). 
[15] O. Biham, I. Furman, V. Pirronello, and G. Vidali, As- 

trophys. J. 553, 595 (2001). 
[16] C. Escudero, Phys. Rev. E 74, 010103(R) (2006). 
[17] M. Assaf and B. Meerson, Phys. Rev. E 74, 041115 

(2006). 

[18] B.C. Gates, JR. Katzer and G.C.A. Schuit, Chemistry 
of Catalytic Processes (McGraw-Hill, New York, 1979). 

[19] SR. Byron, J. Chem. Phys. 30, 1380 (1959); 44, 1378 
(1966). 

[20] For an odd number of particles there is no extinction. 
Here ci = and C2 = 1, and one obtains a true 
steady state: P n = (4{l) n/2 F(n/2) /[nn\ erfi(VQ)], n = 
1,3,5,.... 

[21] G. B. Arfken, Mathematical Methods for Physicists (Aca- 
demic Press, London, 1985). 

[22] V. Elgart and A. Kamenev (private communication). 

[23] CM. Bender and S.A. Orszag, Advanced Mathematical 
Methods for Scientists and Engineers (Springer, New 
York, 1999). 

[24] M. Abramowitz, Handbook of Mathematical Functions 
(National Bureau of Standards, Washington, 1964). 



